Surface energy balance#

In this exercise, we will use pandas, xarray, MetPy, plotly, and rosely to analyze weather station data from a csv file. We will perform data cleaning, manipulation, and visualization to gain insights into the data and explore its characteristics.

Learning objectives:
  • Get familiar with weather station data
  • Quality assessment of measurement time series
  • Data resampling
  • Visualisation with interactive plots
Prerequisites
  • Basic knowledge of Python, Jupyter Notebooks, and data analysis
  • Familiarity with MetPy, Pandas, and Xarray
  • The additional packages rosely and plotly must be installed
  • A csv file containing weather station data (can also be downloaded here)
import numpy as np
import pandas as pd
import math

# Import the plotly library
import plotly.graph_objs as go
from plotly.subplots import make_subplots

from hu_tools import HU_learning_material

#-----------------------------
# Load the learning materials
#-----------------------------
hu = HU_learning_material()

Solar parameters and transmissivity#

#----------------------
# Parameters
#----------------------
ch = 0.0    # High cloud fraction
cm = 0.0    # Middle cloud fraction
cl = 0.0    # Low cloud fraction
lat = 45.0  # Latitude of location
lon = 10.0  # Longitud of location
UTC = 19.5  # UTC at location

#------------------------------------
# Calc parameters and tranmissivity
#------------------------------------
# Create a date vector for which the solar parameters and transmisivity are calculated
timestamp = pd.date_range(start='01-01-2002', end='31-12-2002', freq='D')

# Calculate transmissivity using th HU-Tools
data  = hu.tau(ch, cm, cl, lat, lon, UTC, timestamp.dayofyear, 173, timestamp.size)

# Local elevation angle using th HU-Tools
angle = np.rad2deg(np.arcsin(hu.local_elevation_angle(lat,lon,UTC,timestamp.dayofyear,173,timestamp.size)))

# Create a pandas DataFrame from the calculated arrays and the date vector
df = pd.DataFrame(data, index=timestamp, columns=['tau'])

# Add the local elevation angle to the DataFrame
df['angle'] = angle

#----------------------
# Plot the results
#----------------------
# Creating the plot with two rows and one column. The plots share the same x-axis so that only labels 
# for the lower plots are shown
fig = make_subplots(rows=2, cols=1, shared_xaxes=True)

# Adding the temperature data as a dashed line in the top panel
fig.add_trace(go.Scatter(x=df.index, y=df.tau, line=dict(color='royalblue', dash='solid'), name='Transmissitivty [-]'),
                        row=1, col=1)

# Adding the precipitation data as a dashed line in the bottom panel
fig.add_trace(go.Scatter(x=df.index, y=df.angle, line=dict(color='green', dash='solid'), name='Local elevation angle [º]'),
                         row=2, col=1)

# Adjusting the layout
fig.update_layout(title='Transmissivity and local elevation angle', plot_bgcolor='white', width=800, height=600,
                  yaxis=dict(title='Transmissitivty [-]', showgrid=True, gridcolor='lightgray', gridwidth=1),
                  yaxis2=dict(title='Local elevation angle [º]', showgrid=True, gridcolor='lightgray', gridwidth=1),
                  xaxis=dict(title='', tickformat='%d.%m.%Y', showgrid=True, gridcolor='lightgray', gridwidth=1),
                  xaxis2=dict(title='Date', tickformat='%d.%m.%Y', showgrid=True, gridcolor='lightgray', gridwidth=1))


# Adjusting the axes
fig.update_xaxes(nticks=10, row=1, col=1)
fig.update_yaxes(nticks=5, row=1, col=1)
fig.update_xaxes(nticks=10, row=2, col=1)
fig.update_yaxes(nticks=5, row=2, col=1)

# Showing the plot
fig.show()
Broaden Knowledge & Comprehension
In what period of time is the sun still above the horizon at 45º north and 10º east after 7:30pm? How sensitive is the transmissivity to clouds in the low, middle and high atmosphere? Look at the equation for transmissivity and find out why there are differences. What is a weakness of the parameterisation?

Shortwave radiation and albedo#

albedo = 0.7  # [-] 
SWin = 342    # W/mˆ2

print('Incoming shortwave radiation: {:.2f}'.format(SWin))
print('Outgoing shortwave radiation: {:.2f}'.format(hu.SWout(SWin,albedo)))
Incoming shortwave radiation: 342.00
Outgoing shortwave radiation: 239.40
Broaden Knowledge & Comprehension
Review the 'Global Energy Fluxes' figure on slide 4 of the presentation and identify the outgoing shortwave radiation. Use this information to determine the global albedo.
import numpy as np
import metpy.calc
from metpy.units import units
import math

# The following functions are needed for radiation method Moelg2009
def solpars(lat):
    """ Calculate time correction due to orbital forcing (Becker 2001)
     and solar parameters that vary on daily basis (Mölg et al. 2003)
     0: day angle (rad); 1: in (deg); 2: eccentricity correction factor;
     3: solar declination (rad); 4: in (deg); 5: sunrise hour angle; 6: day length
    """

    timecorr = np.zeros((366, 4))
    solpars = np.zeros((366, 7))

    for j in np.arange(0, 365):
        # Time correction
        x = 0.9856 * (j + 1) - 2.72
        T2 = -7.66 * math.sin(math.radians(x)) - 9.87 * math.sin(
            2 * math.radians(x) + math.radians(24.99) + math.radians(3.83) * math.sin(math.radians(x)))
        timecorr[j, 0] = j + 1  # Julian Day
        timecorr[j, 1] = x
        timecorr[j, 2] = T2  # Time difference between True Local Time (TLT) and Average Local Time (ALT)
        timecorr[j, 3] = T2 * 15 / 60  # Time difference in deg (15°/h)

        # Solar parameters
        tau = 2 * math.pi * (j) / 365
        solpars[j, 0] = tau
        solpars[j, 1] = tau * 180 / math.pi
        solpars[j, 2] = 1.00011 + 0.034221 * math.cos(tau) + 0.00128 * math.sin(tau) + 0.000719 * math.cos(2*tau) + 0.000077 * math.sin(2 * tau)
        solpars[j, 3] = 0.006918 - 0.399912 * math.cos(tau) + 0.070257 * math.sin(tau) - 0.006758 * math.cos(2*tau) + 0.000907 * math.sin(2 * tau) - 0.002697 * math.cos(3 * tau) + 0.00148 * math.sin(3 * tau)
        solpars[j, 4] = solpars[j, 3] * 180 / math.pi
        solpars[j, 5] = math.acos(-math.tan(lat * math.pi / 180) * math.tan(solpars[j, 3])) * 180 / math.pi
        solpars[j, 6] = 2 / 15 * solpars[j, 5]

    # Duplicate line 365 for years with 366 days
    solpars[365, :] = solpars[364, :]
    timecorr[365, :] = timecorr[364, :]

    return solpars, timecorr


def haversine(lat1, lon1, lat2, lon2):
    """ This function calculates the distance between two points given their longitudes and latitudes
     based on the haversine formula. """

    lat1_rad = math.radians(lat1)
    lat2_rad = math.radians(lat2)
    lon1_rad = math.radians(lon1)
    lon2_rad = math.radians(lon2)
    delta_lat = lat2_rad - lat1_rad
    delta_lon = lon2_rad - lon1_rad
    a = ((math.sin(delta_lat / 2)) ** 2 + math.cos(lat1_rad) * math.cos(lat2_rad) * (math.sin(delta_lon / 2)) ** 2) ** 0.5
    d = 2 * 6371000 * math.asin(a)
    return d


def relshad(dem, mask, lats, lons, solh, sdirfn):
    """ This function calculates the topographic shading based on Mölg et al. 2009
     Input:
           dem:    A DEM of the study region, that also includes surrounding terrain
           mask:   A glacier mask
           lats:   The latitudes
           lons:   The longitudes
           solh:   The solar elevation (degrees)
           sdirfn: The illumination direction (in degrees from north)
     Output:
           illu:   Grid illu containing 0 = shaded, 1 = in sun
    """

    z = dem
    illu = dem * 0.0  # create grid that will be filled
    illu[:, :] = np.nan

    # Define maximum radius (of DEM area) in degreees lat/lon
    rmax = ((np.linalg.norm(np.max(lats) - np.min(lats))) ** 2 + (np.linalg.norm(np.max(lons) - np.min(lons))) ** 2) ** 0.5
    nums = abs(int(rmax * len(lats) / (lats[0] - lats[-1])))

    # Calculate direction to sun
    beta = math.radians(90 - sdirfn)
    dy = math.sin(beta) * rmax  # walk into sun direction (y) as far as rmax
    dx = math.cos(beta) * rmax  # walk into sun direction (x) as far as rmax

    # Extract profile to sun from each (glacier) grid point
    for ilat in np.arange(1, len(lats) - 1, 1):
        for ilon in np.arange(1, len(lons) - 1, 1):
            if mask[ilat, ilon] == 1:
                start = (lats[ilat], lons[ilon])
                targ = (start[0] + dy, start[1] + dx)  # find target position

                # Points along profile (lat/lon)
                lat_list = np.linspace(start[0], targ[0], nums)  # equally spread points along profile
                lon_list = np.linspace(start[1], targ[1], nums)  # equally spread points along profile

                # Don't walk outside DEM boundaries
                lat_list_short = lat_list[(lat_list < max(lats)) & (lat_list > min(lats))]
                lon_list_short = lon_list[(lon_list < max(lons)) & (lon_list > min(lons))]

                # Cut to same extent
                if (len(lat_list_short) > len(lon_list_short)):
                    lat_list_short = lat_list_short[0:len(lon_list_short)]
                if (len(lon_list_short) > len(lat_list_short)):
                    lon_list_short = lon_list_short[0:len(lat_list_short)]

                # Find indices (instead of lat/lon) at closets gridpoint
                idy = (ilat, (np.abs(lats - lat_list_short[-1])).argmin())
                idx = (ilon, (np.abs(lons - lon_list_short[-1])).argmin())

                # Points along profile (indices)
                y_list = np.round(np.linspace(idy[0], idy[1], len(lat_list_short)))
                x_list = np.round(np.linspace(idx[0], idx[1], len(lon_list_short)))

                # Calculate ALTITUDE along profile
                zi = z[y_list.astype(int), x_list.astype(int)]

                # Calclulate DISTANCE along profile
                d_list = []
                for j in range(len(lat_list_short)):
                    lat_p = lat_list_short[j]
                    lon_p = lon_list_short[j]
                    dp = haversine(start[0], start[1], lat_p, lon_p)
                    d_list.append(dp)
                distance = np.array(d_list)

                # Topography angle
                Hang = np.degrees(np.arctan((zi[1:len(zi)] - zi[0]) / distance[1:len(distance)]))

                if np.max(Hang) > solh:
                    illu[idy[0], idx[0]] = 0
                else:
                    illu[idy[0], idx[0]] = 1

    return illu


def LUTshad(solpars, timecorr, lat, elvgrid, maskgrid, lats, lons, STEP, TCART):
    """ This function calculates the look-up-table for topographic shading for one year.
     Input:
           solpars:   Solar parameters
           timecorr:  Time correction due to orbital forcing
           lat:       Latitude at AWS
           elvgrid:   DEM
           maksgrid:  Glacier mask
           lats:      Latitudes
           lons:      Longitudes
           STEP:      Time step (s)
           TCART:     Time correction due to difference MLT - TLT
     Output:
           shad1yr:   Look-up-table for topographic shading for 1 year
    """

    hour = np.arange(1, 25, 1)
    shad1yr = np.zeros((int(366 * (3600 / STEP) * 24), len(lats), len(lons)))  # Array (time,lat,lon)
    shad1yr[:, :, :] = np.nan

    # Go through days of year
    for doy in np.arange(0, 366, 1):

        soldec = solpars[doy, 3]  # solar declination (rad)
        eccorr = solpars[doy, 2]  # eccenctricy correction factor
        tcorr = timecorr[doy, 3]  # time correction factor (deg)

        # Go through hours of day
        for hod in np.arange(0, 24, int(STEP / 3600)):

            # calculate solar geometries
            stime = 180 + (15 / 2) - hod * 15 - tcorr + TCART
            sin_h = math.sin(soldec) * math.sin(lat * math.pi / 180) + math.cos(soldec) * math.cos(lat * math.pi / 180) * \
                    math.cos(stime * math.pi / 180)
            cos_sol_azi = (sin_h * math.sin(lat * math.pi / 180) - math.sin(soldec)) / math.cos(math.asin(sin_h)) / \
                    math.cos(lat * math.pi / 180)

            if stime > 0:
                solar_az = math.acos(cos_sol_azi) * 180 / math.pi
            else:
                solar_az = math.acos(cos_sol_azi) * 180 / math.pi * (-1)

            solar_h = math.asin(sin_h) * 180 / math.pi

            sdirfn = 180 - solar_az

            # Calculation (1 = in sun, 0 = shaded, -1 = night)
            if sin_h > 0.01:
                illu = relshad(elvgrid, maskgrid, lats, lons, solar_h, sdirfn)
                shad1yr[round(doy * (3600 / STEP) * 24 + (hod * 3600 / STEP)), maskgrid == 1] = illu[maskgrid == 1]
            else:
                shad1yr[round(doy * (3600 / STEP) * 24 + (hod * 3600 / STEP)), maskgrid == 1] = -1.0

    return shad1yr


def LUTsvf(elvgrid, maskgrid, slopegrid, aspectgrid, lats, lons):
    """ This function calculates the look-up-table for the sky-view-factor for one year.
     Input:
           elvgrid:   DEM
           maksgrid:  Glacier mask
           slopegrid: Slope
           aspectgrid:Aspect
           lats:      Latitudes
           lons:      Longitudes
    """

    slo = np.radians(slopegrid)
    asp = np.radians(aspectgrid)
    res = elvgrid * 0
    count = 0
        
    # Go through all directions (0-360°)
    for azi in np.arange(10, 370, 10):
        # Go through all elevations (0-90°)
        for el in np.arange(2, 90, 2):
            illu = relshad(elvgrid, maskgrid, lats, lons, el, azi)
            a = ((math.cos(np.radians(el)) * np.sin(slo) * np.cos(asp - np.radians(azi))) + (np.sin(np.radians(el)) * np.cos(slo)))
            a[a < 0] = 0
            a[a > 0] = 1
            a[illu == 0] = 0
            res = res + a
            count = count + 1

    vsky = elvgrid * 0
    vsky[:, :] = np.nan
    vsky[maskgrid == 1] = res[maskgrid == 1] / (36 * 44)

    return vsky


def calcRad(solPars, timecorr, doy, hour, lat, tempgrid, pgrid, rhgrid, cldgrid, elvgrid, maskgrid, slopegrid,
            aspectgrid, shad1yr, gridsvf, STEP, TCART):
    """ This function computes the actual calculation of solar Radiation (direct + diffuse)
     including corrections for topographic shading and self-shading, based on Mölg et al. 2009, Iqbal 1983, Hastenrath 1984.
     Input:
           solpars:   Solar parameters
           timecorr:  Time correction due to orbital forcing
           doy:       Day of year
           hour:      Hour of day
           lat:       Latitude at AWSi
           tempgrid:  Air Temperature
           pgrid:     Air Pressure
           rhgrid:    Relative Humidity
           cldgrid:   Cloud fraction
           elvgrid:   DEM
           maksgrid:  Glacier mask
           slopegrid: Slope
           aspectgrid:Aspect
           shad1yr:   LUT topographic shading
           gridsvf:   LUT Sky-view-factor
           STEP:      Time step (s)
           TCART:     Time correction due to difference MLT - TLT
     Output:
           swiasky:   All-sky shortwave radiation
"""

    # Constants
    Sol0 = 1367          # Solar constant (W/m2)
    aesc1 = 0.87764      # Transmissivity due to aerosols at sea level
    aesc2 = 2.4845e-5    # Increase of aerosol transmissivity per meter altitude
    alphss = 0.9         # Aerosol single scattering albedo (Zhao & Li JGR112), unity (zero) -> all particle extinction is due to scattering (absorption)
    dirovc = 0.00        # Direct solar radiation at overcast conditions (as fraction of clear-sky dir. sol. rad, e.g. 10% = 0.1)
    dif1 = 4.6           # Diffuse radiation as percentage of potenial clear-sky GR at cld = 0
    difra = 0.66         # Diffuse radiation constant
    Cf = 0.65            # Constant that governs cloud impact

    soldec = solPars[doy - 1, 3]  # Solar declination (rad)
    eccorr = solPars[doy - 1, 2]  # Cccenctricy correction factor
    tcorr = timecorr[doy - 1, 3]  # Time correction factor (deg)

    # Output files
    swiasky = elvgrid.copy() + np.nan
    swidiff = elvgrid.copy() + np.nan

    # Mixing ratio from RH and Pres
    mixing_interp = metpy.calc.mixing_ratio_from_relative_humidity(pgrid * units.hPa, tempgrid * units.kelvin, rhgrid * units.percent)
    vp_interp = np.array(metpy.calc.vapor_pressure(pgrid * units.hPa, mixing_interp))

    # Solar geometries
    stime = 180 + (STEP / 3600 * 15 / 2) - hour * 15 - tcorr + TCART
    sin_h = math.sin(soldec) * math.sin(lat * math.pi / 180) + math.cos(soldec) * math.cos(
        lat * math.pi / 180) * math.cos(stime * math.pi / 180)
    mopt = 35 * (1224 * (sin_h) ** 2 + 1) ** (-0.5)
    if sin_h < 0:
        mopt = np.nan

    if sin_h > 0.01:  # Calculations are only performed when sun is there

        # Direct & diffuse radiation under clear-sky conditions
        TOAR = Sol0 * eccorr * sin_h
        TAUr = np.exp((-0.09030 * ((pgrid / 1013.25 * mopt) ** 0.84)) * (
                    1.0 + (pgrid / 1013.25 * mopt) - ((pgrid / 1013.25 * mopt) ** 1.01)))
        TAUg = np.exp(-0.0127 * mopt ** 0.26)
        k_aes = aesc2 * elvgrid + aesc1
        k_aes[k_aes > 1.0] = 1.0  # Aerosol factor: cannot be > 1
        TAUa = k_aes ** (mopt)
        TAUaa = 1.0 - (1.0 - alphss) * (1 - pgrid / 1013.25 * mopt + (pgrid / 1013.25 * mopt) ** (1.06)) * (1.0 - TAUa)
        TAUw = 1.0 - 2.4959 * mopt * (46.5 * vp_interp / tempgrid) / ((1.0 + 79.034 * mopt * (46.5 * vp_interp / tempgrid)) **
                                                                      0.6828 + 6.385 * mopt * (46.5 * vp_interp / tempgrid))
        taucs = TAUr * TAUg * TAUa * TAUw

        sdir = Sol0 * eccorr * sin_h * taucs  # Direct solar radiation on horizontal surface, clear-sky
        Dcs = difra * Sol0 * eccorr * sin_h * TAUg * TAUw * TAUaa * (1 - TAUr * TAUa / TAUaa) / (
                    1 - pgrid / 1013.25 * mopt + (pgrid / 1013.25 * mopt) ** (1.02))  # Diffuse solar radiation, clear sky
        grcs = sdir + Dcs  # Potential clear-sky global radiatio∫n

        # Correction for slope and aspect (Iqbal 1983)
        cos_zetap1 = (np.cos(np.radians(slopegrid)) * np.sin(np.radians(lat)) - np.cos(np.radians(lat)) * np.cos(np.radians(180 - aspectgrid)) *
                      np.sin(np.radians(slopegrid))) * np.sin(soldec)
        cos_zetap2 = (np.sin(np.radians(lat)) * np.cos(np.radians(180 - aspectgrid)) * np.sin(np.radians(slopegrid)) +
                      np.cos(np.radians(slopegrid)) * np.cos(math.radians(lat))) * np.cos(soldec) * np.cos(stime * np.pi / 180)
        cos_zetap3 = np.sin(np.radians(180 - aspectgrid)) * np.sin(np.radians(slopegrid)) * np.cos(soldec) * np.sin(stime * np.pi / 180)
        cos_zetap = cos_zetap1 + cos_zetap2 + cos_zetap3

        # Clear-sky direct solar radiation at surface (aspect & slope corrected)
        swidir0 = Sol0 * eccorr * cos_zetap * taucs
        swidir0[cos_zetap < 0.0] = 0.0  # self-shaded cells set to 0
        illu = elvgrid * 0.0
        illu = shad1yr[int(((doy - 1) * (86400 / STEP)) + (hour / (STEP / 3600))), :, :]
        swidir0[illu == 0.0] = 0.0
        sdir[illu == 0.0] = 0.0

        # Correction for cloud fraction
        swidiff[cldgrid > 0.0] = grcs[cldgrid > 0.0] * (((100 - Cf * 100) - dif1) / 100 * cldgrid[cldgrid > 0.0] +
                                    (dif1 / 100)) * gridsvf[cldgrid > 0.0]  # diffuse amount as percentage of direct rad.
        swidiff[cldgrid == 0.0] = Dcs[cldgrid == 0.0] * gridsvf[cldgrid == 0.0]
        swiasky[:, :] = swidir0 * (1 - (1 - dirovc) * cldgrid) + swidiff  # all-sky solar radiation at surface

    else:
        TOAR = 0.0
        swiasky[maskgrid == 1] = 0 * elvgrid[maskgrid == 1]
        illu = 0.0 * elvgrid - 1

    swiasky_ud = swiasky[::-1, :]
 
    return swiasky_ud
import rasterio
from rasterio.plot import show
from osgeo import gdal
import xarray as xr
from xrspatial import slope, aspect
import rioxarray

dem = xr.open_dataset("./data/paine_dem.tif", engine="rasterio")

dem = dem.sel(y=slice(-50.93,-50.96), x=slice(-73.06,-73.00))
dem.to_netcdf('paine_dem.nc', engine='netcdf4')
#dem['band_data'][0,:,:].plot()

elvgrid = dem['band_data'][0,:,:].data
slopegrid = slope(dem['band_data'][0,:,:]).data
aspectgrid = aspect(dem['band_data'][0,:,:]).data


maskgrid = np.ones_like(aspectgrid)
lons = dem['band_data'].x.data
lats = dem['band_data'].y.data


#solpars_param, timecorr = solpars(lats[0])

# Sky view factor
#svf = LUTsvf(np.flipud(elvgrid), np.flipud(maskgrid), np.flipud(slopegrid), np.flipud(aspectgrid), lats[::-1], lons)
#print('Look-up-table sky view factor done')

# Topographic shading
#shad1yr = LUTshad(solpars_param, timecorr, lats[0], np.flipud(elvgrid), np.flipud(maskgrid), lats[::-1], lons, 3*3600, 0.0)
#print('Look-up-table topographic shading done')
import netCDF4 as nc

dtstep = 3*3600
tcart = 0

# Save look-up tables
Nt = int(366 * (3600 / dtstep) * 24)  # number of time steps
Ny = len(lats)  # number of latitudes
Nx = len(lons)  # number of longitudes

f = nc.Dataset('./data/LUT_Rad.nc', 'w')
f.createDimension('time', Nt)
f.createDimension('lat', Ny)
f.createDimension('lon', Nx)

LATS = f.createVariable('lat', 'f4', ('lat',))
LATS.units = 'degree'
LONS = f.createVariable('lon', 'f4', ('lon',))
LONS.units = 'degree'

LATS[:] = lats
LONS[:] = lons

shad = f.createVariable('SHADING', float, ('time', 'lat', 'lon'))
shad.long_name = 'Topographic shading'
shad[:] = shad1yr

SVF = f.createVariable('SVF', float, ('lat', 'lon'))
SVF.long_name = 'sky view factor'
SVF[:] = svf

f.close()
dem = xr.open_dataset("./data/paine_dem.tif", engine="rasterio")
dem = dem.sel(y=slice(-50.93,-50.96), x=slice(-73.06,-73.00))

dem.to_netcdf('./data/paine_dem.nc')

elvgrid = dem['band_data'][0,:,:]
slopegrid = slope(dem['band_data'][0,:,:])
aspectgrid = aspect(dem['band_data'][0,:,:])


maskgrid = xr.zeros_like(aspectgrid)
lons = elvgrid.x
lats = elvgrid.y

solpars_param, timecorr = solpars(lats[0])

# Calculate solar Parameters
solPars, timeCorr = solpars(lats[0])

print('Read in look-up-tables')
ds_LUT = xr.open_dataset('./data/LUT_Rad.nc')
shad1yr = ds_LUT.SHADING.values
svf = ds_LUT.SVF.values
Read in look-up-tables
import plotly.graph_objects as go
import plotly.subplots as sp
import numpy as np
import xarray as xr
from xrspatial import slope, aspect

#------------------------------
# Read digital elevation model
#------------------------------
dem = xr.open_dataset('./data/paine_dem.nc')

# Read the elevation data from netcdf file
elvgrid = dem['band_data'][0,:,:]
# Read slope data
slopegrid = slope(dem['band_data'][0,:,:]).data
# Read aspect data
aspectgrid = aspect(dem['band_data'][0,:,:]).data

print('Read in look-up-tables')
ds_LUT = xr.open_dataset('./data/LUT_Rad.nc')
shad1yr = ds_LUT.SHADING.values
svf = ds_LUT.SVF.values

# Create additional arrays which are needed by the radiation module
# Temperature grid
tempgrid = np.ones_like(elvgrid) * 285.0
# Pressure grid
pgrid = np.ones_like(elvgrid) * 1000.0
# Relative humditiy grid
rhgrid = np.ones_like(elvgrid) * 70
# Cloud grid
ngrid = np.ones_like(elvgrid) * 0.0
# Mask grid (in case only part of the DEM is used
maskgrid = xr.zeros_like(elvgrid)

# Timestamp
lat = -50.9 # Latitude of the domain center
doy = 15    # Doy of the year
hour = 13   # Hour of the day
tcart = 0   
dtstep = 3*3600

# Get lon/lat values from the elevation model
lons = elvgrid.x
lats = elvgrid.y

#---------------------------
# Calculate solar Parameters
#---------------------------
solPars, timeCorr = hu.solpars(lats[0])

#---------------------------
# Calc radiation
#---------------------------
hu
A = hu.calcRad(solPars, timeCorr, doy, hour, lat, tempgrid[::-1, :], pgrid[::-1, :], rhgrid[::-1, :], 
        ngrid[::-1, :], np.flipud(elvgrid), np.flipud(maskgrid), np.flipud(slopegrid), np.flipud(aspectgrid), shad1yr, svf, dtstep, tcart)

#---------------------------
# Plot data
#---------------------------
# Generate sample data
x = lons
y = lats
X, Y = np.meshgrid(x, y)

# Create the subplots
fig = sp.make_subplots(rows=1, cols=2)

# Add the first subplot
fig.add_trace(go.Contour(x=x, y=y, z=elvgrid, colorscale='Viridis', showscale=True,
                         colorbar=dict(x=0.2, y=-0.22, len=0.4, yanchor='bottom', orientation='h',
                                      title=r'$\text{Incoming shortwave radiation}~\text{W} m^{-2}$', titleside='top')), row=1, col=1)

# Add the second subplot
fig.add_trace(go.Contour(x=x, y=y, z=A, colorscale='YlOrRd', showscale=True,
                         colorbar=dict(x=0.8, y=-0.22, len=0.4, yanchor='bottom', orientation='h', 
                                       title=r'$\text{Incoming shortwave radiation}~\text{W} m^{-2}$', titleside='top')), row=1, col=2)



# Add title and axis labels
fig.update_layout(
    xaxis=dict(title='Longitude',title_standoff=5),
    xaxis2=dict(title='Longitude',title_standoff=5),
    yaxis_title='Latitude',
    width=1200,  # Set the width of the figure
    height=600  # Set the height of the figure
)

# Show the plot
fig.show()
Read in look-up-tables

Longwave radiation#

# Define a emissivity
epsilon = 1.0

# Radiated energy at zero degrees celcius
print('Radiated energy at zero degrees celcius: {:.4}'.format(hu.LW(epsilon, 273)))

# Wave length of emission maximum
print('Wave length of emission maximum: {:.4}'.format(hu.lambda_max(273)))
Evaluate and Analyze
  • Calculate the long-wave radiation for different temperatures. What is the relationship between temperature and radiated energy? [Hint: Take a look at the Stephan-Boltzmann Law]
  • Plot the wavelength of the emission maximum against the radiated energy in the range from 250 K to 310 K. Which gases absorb in this spectral range? Comment on the influence of the most important gases on the greenhouse effect?
Broaden Knowledge & Comprehension
What must be taken into account when measuring long-wave radiation in order to obtain reliable values?

Ground heat flux and heat equation#

#-----------------------
# Boundary conditions
#-----------------------
Ts = 20.   # Temperature at the surface [K]
Tb = 5.    # Temperature at the bottom [K]
D  = 2.    # Depth [m]
Nz = 100   # Number of grid points in the vertical (soil)
dt = 60.   # time step [s]
t7200  = 7200  # integration time [s]
t86400 = 86400 # integration time [s]
alpha = 1.2e-6 # soil thermal diffusivity [m^2 s^-1]

#-----------------------
# Run the heat equation
#-----------------------
T, dz = hu.heat_equation_indices(Ts, Tb, D, Nz, t7200, dt, alpha)
T2, dz = hu.heat_equation(Ts, Tb, D, Nz, t86400, dt, alpha)

#-----------------------
# Make plots
#-----------------------
# Creating the plot with two rows and one column. The plots share the same x-axis so that only labels 
# for the lower plots are shown
fig = make_subplots(rows=1, cols=1, shared_yaxes=True)

# Adding the temperature data as a dashed line in the top panel
fig.add_trace(go.Scatter(x=T, y=-dz*np.arange(Nz), line=dict(color='royalblue', dash='solid'), name='Temperature after 7200 s'),
                        row=1, col=1)

# Adding the temperature data as a dashed line in the top panel
fig.add_trace(go.Scatter(x=T2, y=-dz*np.arange(Nz), line=dict(color='green', dash='solid'), name='Temperature after 86400'),
                        row=1, col=1)

# Adjusting the layout
fig.update_layout(title='Soil temperature', plot_bgcolor='white', width=800, height=600,
                  yaxis=dict(title='Depth [m]', showgrid=True, gridcolor='lightgray', gridwidth=1),
                  xaxis=dict(title='Temperature [K]', tickformat='%d.%m.%Y', showgrid=True, gridcolor='lightgray', gridwidth=1))


# Adjusting the axes
fig.update_xaxes(nticks=10, row=1, col=1, range=[-2,20], dtick=5)
fig.update_yaxes(nticks=5, row=1, col=1)

# Showing the plot
fig.show()
Evaluate and Analyze:
  • What happens if the integration time is set to a very long period?
  • How can the model be improved? What influence does soil thermal diffusivity have on the temperature profile?
  • Take diffusivity values for different soil types from the literature and compare the different temperature profiles.
Broaden Knowledge & Comprehension
How can the soil heat flux be determined from the temperature profile?
#-----------------------------
# Create pandas dataframe with
# the surface temperature
#-----------------------------
# timestamp from/to
timestamp = pd.date_range(start='2001-01-01T00:00:00', end='2003-12-31T00:00:00', freq='D')

# Surrogate surface temperature timeseries
T = 10 - 20 * np.sin((2*math.pi*timestamp.dayofyear)/365)

# Create a pandas DataFrame from the calculated arrays
df = pd.DataFrame(T, index=timestamp, columns=['T'])

#-----------------------------
# Run the time-dependent heat
# equation
#-----------------------------
D  = 20   # Depth of the domain
dz = 0.5  # Spacing between grid points
T, D, dz = hu.heat_equation_daily(df['T'], D, dz, alpha=1.2e-6)

#-----------------------------
# Plot the results
#-----------------------------
# Generate the plotting grid
# y-axis values
y_values = np.arange(-D,0,dz)
# then the x-axis values
x_values = np.arange(len(df.index))
# Generate the 2D plotting grid
X, Y = np.meshgrid(x_values, y_values)

# Create the filled contour trace
trace = go.Contour(
    x=df.index,
    y=y_values,
    z=T[::-1],
    colorscale='Viridis',
    contours=dict(
        coloring='fill',
        showlabels=True
    ),
    colorbar=dict(
        title=dict(
            text='Temperature [ºC]',
            #standoff=15,
            side='right'
        )
    )
)

# Create the layout
layout = go.Layout(
    title='Soil temperature',
    xaxis=dict(title='Date', tickformat='%d.%m.%Y', tickangle=45, dtick='M3'),
    yaxis=dict(title='Depth [m]'),
    width=1000, height=500,
)

# Create the figure
fig = go.Figure(data=[trace], layout=layout)

# Show the plot
fig.show()
Broaden Knowledge & Comprehension
Change the soil thermal diffusivity and see how the profile changes.
# Creating the plot with two rows and one column. The plots share the same x-axis so that only labels 
# for the lower plots are shown
fig = make_subplots(rows=1, cols=1, shared_xaxes=True)

# Adding the temperature data as a dashed line in the top panel
fig.add_trace(go.Scatter(x=df.index, y=T[10,:], line=dict(color='royalblue', dash='solid'), name='5 m'),
                        row=1, col=1)

# Adding the temperature data as a dashed line in the top panel
fig.add_trace(go.Scatter(x=df.index, y=T[20,:], line=dict(color='orange', dash='solid'), name='10 m'),
                        row=1, col=1)

# Adjusting the layout
fig.update_layout(title='Soil temperature', plot_bgcolor='white', width=1000, height=600,
                  yaxis=dict(title='Soil temperature [K]', showgrid=True, gridcolor='lightgray', gridwidth=1),
                  xaxis=dict(title='', tickformat='%d.%m.%Y', showgrid=True, gridcolor='lightgray', gridwidth=1))


# Adjusting the axes
fig.update_xaxes(nticks=10, row=1, col=1, tickformat='%d.%m.%Y', tickangle=45, dtick='M3')
fig.update_yaxes(nticks=5, row=1, col=1)

# Showing the plot
fig.show()
help(hu.emissivity)

Load weather station data#

In this example, the code reads in the 30-min weather stations data from a CSV file using the ‘read_csv()’ method from pandas, with the parse_dates and index_col parameters set to True and 0, respectively. This ensures that the first column of the CSV file (the timestamps) is used as the index of the DataFrame, and that the timestamps are converted to datetime objects.

import pandas as pd
import numpy as np

# Load CSV file
df = pd.read_csv("https://raw.githubusercontent.com/sauterto/clim_env_hydro/main/docs/nb/data/FLX_CH-Dav.csv", parse_dates=True, index_col=0)

# The file contains:
# Air temperature       :: t2m
# Relative humdity      :: RH
# Precipitation         :: precip
# Wind speed            :: WS
# Wind direction        :: WD
# Net radiation         :: NETRAD
# Incoming shortwave    :: SW_IN
# Outgoing shortwave    :: SW_OUT
# Incoming longwave     :: LW_IN
# Outgoing longwave     :: LW_OUT
# Sensible heat flux    :: H
# Latent heat flux      :: LE
# Ground heat flux      :: QG
# Surface temperature   :: TS
Broaden Knowledge & Comprehension
Which other file formats can be read in with pandas?

Replace missing data#

Next, we want to view the data to ensure it has been loaded correctly. We use the head() method from Pandas to display the first five rows of the DataFrame.

The DataFrame has missing values with the dummy value -9999. We then use the replace method to replace all occurrences of -9999 with NaN (Not a Number), using the NumPy np.nan constant. Finally, we print the updated DataFrame. The inplace=True argument ensures that the original DataFrame is modified, rather than a copy being created.

# replace -9999 with NaN
df.replace(-9999, np.nan, inplace=True)

# Keep only the surface temeperature
ds = df['TS']

# check for missing values
ds.dropna(inplace=True)

# Show time series
ds
Hint! Compare the SW_OUT and LW_OUT column with those of the previous table.

To check for missing values in a Pandas DataFrame, you can use the isna() or isnull() method, which returns a Boolean DataFrame of the same shape indicating which cells contain missing values.

You can also use the isna() or isnull() method along with the sum() method to count the number of missing values in each column:

# resample dataframe to monthly mean values with different aggregation for each column and return NaN for 
# columns with insufficient valid elements
threshold =  50 # What is the maximum percentage of NaNs that may be included in the averaging? Here, we set the threshold to 10%

# This is a pythonic way in solving this problem
ds_daily = ds.resample('1D').agg({'TS': lambda x: x.dropna().mean() if (((x.isna().sum())/len(x))*100) < threshold else np.nan})

# Find dates with missing data
missing_dates = ds_daily['TS'][ds_daily['TS'].isnull()].index
print("Dates with missing data:\n",ds_daily.loc[missing_dates])

# Check if the time series is continuous
if pd.date_range(start=ds_daily.index.min(), end=ds_daily.index.max(), freq='D').difference(ds_daily.index).empty:
    print('The time series is continuous \n')
else:
    print('The time series is not continuous \n')
    
# Plot the temperature time series
ds_daily['TS'].plot()
#-----------------------------
# Load the learning materials
#-----------------------------
hu = HU_learning_material()

#-----------------------------
# Run the time-dependent heat
# equation
#-----------------------------
D  = 10   # Depth of the domain
dz = 0.5  # Spacing between grid points

df_hourly_subset = ds_daily.loc['2010':'2014']

T, D, dz = hu.heat_equation_daily(df_hourly_subset['TS'], D, dz, alpha=1.2e-6)

#-----------------------------
# Plot the results
#-----------------------------
# Generate the plotting grid
# y-axis values
y_values = np.arange(-D,0,dz)
# then the x-axis values
x_values = np.arange(len(df_hourly_subset.index))
# Generate the 2D plotting grid
X, Y = np.meshgrid(x_values, y_values)

# Create the filled contour trace
trace = go.Contour(
    x=df_hourly_subset.index,
    y=y_values,
    z=T[::-1],
    colorscale='Viridis',
    contours=dict(
        coloring='fill',
        showlabels=True
    ),
    colorbar=dict(
        title=dict(
            text='Temperature [ºC]',
            #standoff=15,
            side='right'
        )
    )
)

# Create the layout
layout = go.Layout(
    title='Soil temperature',
    xaxis=dict(title='Date', tickformat='%d.%m.%Y', tickangle=45, dtick='M3'),
    yaxis=dict(title='Depth [m]'),
    width=1000, height=500,
)

# Create the figure
fig = go.Figure(data=[trace], layout=layout)

# Show the plot
fig.show()